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We introduce an approach to derive realistic Coulomb interaction terms in free standing layered 
materials and vertical heterostructures from ab-initio modelling of the corresponding bulk materials. 

To this end, we establish a combination of calculations within the framework of the constrained 
random phase approximation, Wannier function representation of Coulomb matrix elements within 
some low energy Hilbert space and continuum medium electrostatics, which we call Wannier function 
continuum electrostatics (WFCE). For monolayer and bilayer graphene we reproduce full ab-initio 
calculations of the Coulomb matrix elements within an accuracy of 0.2 eV or better. We show that 
realistic Coulomb interactions in bilayer graphene can be manipulated on the eV scale by different 
dielectric and metallic environments. A comparison to electronic phase diagrams derived in [M. M. 

Scherer et al., Phys. Rev. B 85, 235408 (2012)] suggests that the electronic ground state of bilayer 
graphene is a layered antiferromagnet and remains surprisingly unaffected by these strong changes 
in the Coulomb interaction. 


I. INTRODUCTION 

Since the isolation of monolayer graphene^, the library 
of experimentally available two dimensional (2D) mate¬ 
rials has been continuously growing^ and includes now 
e.g. graphene, hexagonal boron nitride, silicene or the 
class of transition metal dichalcogenides (TDMCs). To¬ 
gether with recent developments in on-demand stacking 
of these atomically thin crystals, a whole new class of 
hybrid materials is coming into reach, nowA While these 
2D structures hold promises from the realization of exotic 
electronic quantum phases at elevated temperatures^ to 
concepts for optoelectronic information processing and 
light harvesting^, their theoretical description is very 
challenging as they generally combine structural com¬ 
plexity with pronounced electronic interaction effects 
such as renormalized quasiparticles or competing elec¬ 
tronic phases^. 

One general theoretical strategy towards realistic de¬ 
scriptions of interacting electron systems is the combi¬ 
nation of ab-initio and model Hamiltonian approaches, 
where single particle and interaction parameters are de¬ 
rived from first principles calculations^. Here, one typ¬ 
ically considers quantum lattice models like extended 
multiband Hubbard models, which describe electrons 
within some set of low energy bands interacting with each 
other. Realistic Coulomb interaction matrix elements en¬ 
tering these models should be appropriately screened, i.e. 
account for screening due to those states which are not 
exciplicitly treated in the low energy models, and can be 
calculated from first principles using the so-called con¬ 
strained Random Phase Approximation (cRPA)ii. The 
computational demand of these calculations is compara¬ 
ble to GW calculations, which makes the treatment of 
complex heterostructures e.g. in plane wave based ap¬ 


proaches very challenging. 

Essentially two different strategies have been devel¬ 
oped to circumvent computational problems in obtain¬ 
ing appropriately screened interactions for thin films or 
layered materials. First, in long-wavelength approaches 
to layered materials model dielectric functions based on 
a description of the two-dimensional screening in terms 
of macroscopic electrodynamics can be straightforwardly 
employed (see e.g. Refs. EH and EH). Second, modified 
Coulomb interactions involving e.g. a truncation in the 
vertical direction^ 3 or unscreening in terms of model di¬ 
electric functions^ can be employed directly on the fully 
microscopic GW level to reach faster convergence of the 
screened interactions in repeated slab approaches. While 
the first set of approaches comes at the advantage of al¬ 
most no computation cost in obtaining screened Coulomb 
interactions, it generally relies on a-priori unknown ad¬ 
justable parameters. The second set of approaches con¬ 
tains all microscopic real material informations but re¬ 
quires a new fully microscopic calculation when the di¬ 
electric environment of some layered material (e.g. the 
substrate) is changed and remains still computationally 
demanding. 

In this paper, we introduce a bridge between these two 
complementary classes of approaches and develop an ap¬ 
proximate very simple yet accurate approach to derive 
realistic Coulomb interaction matrix elements for elec¬ 
trons in free standing layered materials and vertical het¬ 
erostructures from cRPA modelling of the correspond¬ 
ing bulk materials. To this end, we combine Wannier 
function representations of the Coulomb matrix elements 
within some low energy Hilbert space of interest with con¬ 
tinuum medium electrostatics, as we explain in section 
IIII1 This allows us to avoid repeated slab calculations on 
the cRPA level. In section HVl we illustrate our Wannier 
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function based approach with the example of graphene, 
bilayer graphene as well as related heterostructures like 
Ir intercalated graphene. A particular advantage of the 
approach introduced, here, is that one can very easily 
assess how different environments affect Coulomb inter¬ 
actions in realistic layered materials, as we show with the 
example of bilayer graphene in section [V] 

II. MODEL HAMILTONIAN AND THE CRPA 
APPROACH 

For the graphene based systems to be studied, here, we 
consider an effective model, which includes the carbon 
p z orbitals (i.e. the n bands) and treats the a bands as 
well as states at higher energies as the ’’rest”—. We can 
describe such a system with a generalized Hubbard model 
for the p z orbitals with the many-body Hamiltonian 

H = — t Y + lJ 00 Y n n n il (!) 

<i,j>,cr i 

+ — ^ ' U i 

i#j , 

cr,cr 

where c- la annihilates an electron with spin a £ {f, ),} 
at site i and = c? c^. The index i = (i, A or B) 
labels the sublattice (A, B) and the unit cell centered at 
position Rj. 

The cRPAii is used to derive the effective partially 
screened Coulomb interaction terms Uij entering the 
model of Eq. CD from first principles. Within the stan¬ 
dard random phase approximation (RPA) the screened 
Coulomb interaction is defined as W{ q, w) = [1 — 
w(q)P(q, u;)] _1 ?;(q) witla the bare Coulomb repulsion 
i>(q) and the momentum and frequency dependent po¬ 
larization function P(q, ui). We note that all quantities 
are non-local in space, e.g., v(r, r', q), where r and r' are 
restricted to the unit cell. Using a suitable basis the func¬ 
tions become matrices and the equation above becomes 
a matrix equation with 1 being the unit matrix. 

To obtain the effective Coulomb interaction, the func¬ 
tion P(q, uj) is divided into a P n and a P r part. The 
former accounts for charge fluctuations that are induced 
by virtual transitions within the 7r bands only, i.e., from 
7T to 7r*. All other transitions are included in the re¬ 
mainder P r . The latter give rise to screening effects that 
reduce the effective Coulomb interaction between the 7r 
electrons of the model 

U(q,u) = - --- = e" 1 (q,w)u(q), (2) 

1 - u(q)P r (q,w) 

where e(q, w) is the matrix representation of the micro¬ 
scopic dielectric function e(r, r', q, w) and e -1 its inverse. 
The interaction U is in general momentum and frequency 
dependent and can be long ranged. 

We evaluate these bare and screened Coulomb ma¬ 
trix elements in an all-electron mixed product basi s 17 ’ 18 


based on the full-potential linearized augmented-plane- 
wave (FLAPW) methodic—. For this purpose we start 
with a density functional theory (DFT) calculation em¬ 
ploying the FLEUR cod e 22 ’ 23 to obtain the corresponding 
ground states within the FLAPW method. Afterwards, 
we use the SPEX cod o 17 ’ 18 to calculate the bare and 
screened Coulomb matrix elements in the constrained 
random phase approximation 2 ^. Technical details con¬ 
cerning these calculations can be found in the appendix 
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III. COMBINATION OF WANNIER 
FUNCTIONS AND CONTINUUM 
ELECTROSTATICS 

In this section, we explain our approach to derive 
appropriately screened Coulomb interaction matrix el¬ 
ements for electrons in two dimensional materials (e.g. 
graphene) and their heterostructures on the basis of 
Coulomb interaction matrix elements from parent three 
dimensional bulk systems (e.g. graphite). To this end, 
we recall first the continuum electrodynamic description 
of layered materials, free standing monolayers as well as 
heterostructures. Afterwards, the continuum formula¬ 
tion is embedded into a quantum lattice description in 
terms of localized Wannier functions. 


A. Continuum electrostatic description 



Figure 1. (Color online) Left: Three dimensional layered 
material. The infinite stack of single layers is numbered by j. 
Right: Corresponding two dimensional monolayer in a vari¬ 
able dielectric surrounding. The arising electric field of two 
charges q a and qb is indicated by dashed lines. The dielectric 
properties in the layer is denoted by ei, which is actually a 
non-local and thus qy dependent function (see main text). 


The problem of screening in terms of continuum elec¬ 
trostatics in a layered material is illustrated in Fig. |T| 
and has been considered in Refs. l25l - |27i . We generalize 
these works to include non-local screening effects. To this 
end, we start with the Poisson equation in the presence 
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of matter: 

V [ dV [<5(r - r') + X(r - r')] E(r') = —p ext (r), (3) 

J £o 

which relates the divergence of the dielectric displace¬ 
ment (left) to the density of the external, i.e., unbound, 
charge (right). The dielectric displacement is written in 
terms of the electric field E(r') and the induced polar¬ 
ization f X(r — r')E(r , )d 3 r', where the susceptibility X 
is generally non-local but depends only on differences of 
the spatial coordinates r — r'. 

These coordinates can be separated into the in-plane 
and vertical components according to r = (r || , z): 

V yd 3 r' [<5(r|| — r',, )S(z — z') + X(r,, — r'||, z — z')\ (4) 

■ E(r'||,z') = — Pext(r\\,z). 

£o 

The electric field E(r'||, z') can be obtained from the elec¬ 
trostatic potential E(r'||,z') = — V$(r'||, z'). 

Due to the translational invariance of the continuum 
medium, Eq. (j4j) becomes diagonal in Fourier space and 
shows a simple algebraic relation between the suscepti¬ 
bility, the potential, and the external charge density 

q 2 [i + x(q||,g*)] $(q|i,&) = —Pext(q||,9«)- (5) 

£o 

Instead of the susceptibility we can consider the dielectric 
function e(qj|, q z ) = 1+X(q||, q z ) and express the solution 
of the Poisson equation © in the form 


function e(r, r', qy, q z ) so that we may assume e(q||,g-) 
to be a scalar function. The embedding is not as easily 
done as in Refs. [ 25 H 27 L because we face the problem of 
non-localities, in particular the q z dependence, which de¬ 
scribes the periodicity of the bulk material in z-direction. 
Clearly, an assumption on how non-localities translate 
from bulk to monolayer or heterostructures has to be 
made. Here, we continue with the simplest possible ap¬ 
proximation and neglect all non-localities in z-direction, 
i.e. we replace 

h r ,h 

e(q||,&) -> £i(q||) = TT / dg z e(q||,g z ), (8) 
Z7T J-Tt/h 

where h plays the role of an effective layer thickness. This 
definition is plausible as the two-dimensional embedding 
breaks the periodicity in z-direction, and only the local 
term of £ should remain relevant. In this sense, Eq. © 
gives this local term as the Fourier transformation to the 
center of the monolayer, i.e., z = 0. A similar formula 
was used in Ref. 0 to define a ’’two-dimensional macro¬ 
scopic dielectric function”. 

We now assume that £i(qj|) is constant on the whole 
width of the monolayer, i.e., for |z| < h/ 2 , and consider 
two dielectric materials on both sides with dielectric con¬ 
stants e 2 and £3 according to Fig. Q] which gives 

f £2 z > f 

e (qii>2) = s £ i(qii) M < I ■ ( 9 ) 

1 £3 z < =r 


, s 1 . , 3>ext(q||, <7z) 

(q|l ’ 9z) - £oe(q||,gz)<7 2/5ext(q||,9z) - £( q|l , 9z ) ’ 

( 6 ) 

where we introduced 3> ex t(q||, <7z), which is the electro¬ 
static potential created by the external charge. 

In turn, Eq. © can be regarded as a definition of the 
dielectric function: 


£ (q|b9z) 


4 > ext(q||; qz) 

$(q||»9*) 


( 7 ) 


On the other hand, the dielectric function can be 
obtained from first principles, for example using 
the random-phase approximation. In this context, 
^ext(q| ■ qz) plays the role of the bare Coulomb poten¬ 
tial, which is responsible for the bare Coulomb inter¬ 
action u 3D (q||, q z ) = —e4> ext (q||, q~) and $(qy ,q z ) com¬ 
prises u 3D (q||,g z ) and the potential created by the in¬ 
duced charge. So, $(q||,< 7 z) corresponds to the screened 
interaction C7 3D (q||, q z )- 

We now proceed to relate the bulk dielectric function 
to the dielectric function of a two-dimensional system 
embedded into a dielectric environment, as shown in Fig. 
|TJ We assume that the former has been determined from 
first principles. As we will explicate later, we only mod¬ 
ify the leading eigenvalue of the microscopic dielectric 


To find the appropriately screened interaction [/ 2D (q||) 
between two electrons in the central layer, we consider 
the electrostatic problem of an oscillating two dimen¬ 
sional charge density p(q||, z) = p 2D (q\\)S(z) in the center 
of the monolayer, which is at z = 0. The resulting elec¬ 
trostatic potential $(q||, z) is not the same as in the bulk 
[Eq. ©] due to modified screening and the fact that 
we are now considering a 2D charge density confined to 
z = 0. 

In the dielectric continuum model, the modification 
of the screening is caused by the formation of image 
charges at the interfaces between the dielectrics. Let 
us first consider a single interface, where the dielectric 
constant changes from t\ to £ 2 , and a test charge q in 
region 1. For an observer in region 1, the induced po¬ 
larization has the form of an image charge of magnitude 
<?(ei — £ 2 )/(ei + £ 2 ) that is located in region 2 at an equal 
distance from the interface as the test charge^. If one 
has more than one interface, as in the case of our dielec¬ 
tric model, the charge is reflected infinitely many times 
(as light between two parallel optical mirrors), giving rise 
to an infinite number of image charges of ever decreasing 
magnitude perpendicular to the interfaces. In the more 
general case of the two-dimensional charge distribution, 
each Fourier component of p(r) is mirrored according to 
the corresponding component of the dielectric function 
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and the ratio is 




gi(qn) ~ £ j 

£l(q||) + £j 


with j = 2 and j = 3 for the upper and lower inter¬ 
face, respectively. Several works have treated this sit¬ 
uation, a dielectric monolayer sandwiched between two 
semi-infinite dielectric materials^ - —. We use here a spe- 
(10) cial case of a formula derived in Ref. [23, giving the ef¬ 
fective two-dimensional dielectric function 

I 


£eff(q||) 


_£i(q|j) [l -£ 2 (q||)£3(q||)e 2<?l|/l ]_ 

1+ [£ 2 (q||) +£ 3 (q||)] e"«n ft '-|-£2(q||)£3(q||)e“ 2,? ii /t ' 


( 11 ) 


With this equation, the two-dimensional screened inter¬ 
action (in the long-wavelength limit, see below) can be 
written as 


^ 2 D (qn) 


^ 2D (qn) 

£ eff(qn)’ 


( 12 ) 


where f 2 D (q|[) = i' 3 D (q||,z = 0 ) is the bare interaction 
in the 2D system. 

By evaluating Eq. m with the help of Eq. © we are 
now able to calculate the screened interaction between 
electrons in the free standing or embedded two dimen¬ 
sional monolayer directly from the three dimensional lay¬ 
ered bulk properties with the main approximation being 
the neglect of all non-localities of dielectric response in 
the vertical direction. We will assess the quality of this 
approximation in section IIV1 


B. Wannier function based formulation 

The generalized Hubbard model, Eq. Q, involves ma¬ 
trix elements in terms of Wannier functions^ Thus, the 
continuum electrostatic description developed in the pre¬ 
vious section has to be transferred to this Wannier func¬ 
tion basis. 

In a situation with one Wannier orbital per unit cell 
and density-density type interactions only, the interac¬ 
tion terms of Eq. © can be expressed in fc-space, 
f7(q)n q n_ q , with the Fourier transformed electron 

density ?i q = ]Tk a e* qRi nicr and the matrix elements 
U{ q) = e' qR * t/o.i ^ This kind of transformation can 
be made for the bare i> 3 D (q) and the partially screened 
matrix elements U( q). It is then natural to identify 
these matrix elements with the bare and screened in¬ 
teractions discussed in the previous section and to de¬ 
rive the screened Coulomb interaction of a free standing 
monolayer from its three dimensional layered bulk coun¬ 
terpart simply according to the algorithm summarized in 
Eqs. 0 to (fill) . 

Often, model descriptions of real materials like 
graphene or transition metal dichalcogenides are more 
complicated as they involve multiple Wannier orbitals per 
unit cell and also general non-density-density interaction 


terms. Then we have to deal with full matrix representa¬ 
tions of the interactions. For example, U( q) is a matrix 
with the elements g ( q), which depend in general on 

two initial momenta k, k' and the momentum transfer q 
as well as four orbital indices a, (3, 7 and S. 

In this case further approximations are helpful to de¬ 
rive from the bulk Coulomb interaction the correspond¬ 
ing monolayer terms. First, we neglect Coulomb assisted 
hopping terms between sites, which means tracing out 
the k and k' dependencies. Then, the Coulomb inter¬ 
action depends on momentum transfer q but not on the 
initial momenta k and k'. 

To combine the macroscopic electrostatic description 
of the previous section with the representation in a Wan¬ 
nier basis, we still have to account for the orbital de¬ 
pendencies. We can represent the bare and screened 
Coulomb interaction as well as the dielectric function as 
quadratic matrices, using generalized indices a = {a , d} 
and /3 = {/3,7}: U a p 7 s( q) = tTjfq), which can be inter¬ 
preted as interaction energies of generalized charge den¬ 
sity waves n & (q) = )C k (cjj(k+q)c, 5 (k)) and rig(-q). The 
leading eigenvalue and eigenvector of U & p( q) corresponds 
to the charge density modulations with the longest wave¬ 
length, i.e. those charge density waves where screen¬ 
ing effects due to the environment are supposed to be 
strongest. By definition, continuum medium electrostat¬ 
ics describes the long wavelength response of a medium 
to long wavelength electric field / potential variations. 
We thus correct the leading eigenvalue of U-s( q) accord¬ 
ing to the algorithm from the previous section, while we 
assume for all other eigenvalues the same screening as 
in the bulk. The full algorithm which we refer to as 
“Wannier Function Continuum Electrostatic” (WFCE) 
approach can be divided into several steps, which are il¬ 
lustrated in the flowchart of Fig. [2] 

We start with the bare, v? R (q), and partially screened 
Coulomb interactions f7??(q) obtained from the ab-initio 

otp v ' 

calculations for the three dimensional bulk of the layered 
material. We then diagonalize u 3 R (q), obtain the corre¬ 
sponding transformation matrix S’(q), and assume that 
S(q) also diagonalizes C/??(q). In this way, we define the 
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Figure 2. Flowchart of the Wannier function continuum electrostatics (WFCE) algorithm to obtain the screened Coulomb 
matrix elements of a freestanding monolayer directly from the Coulomb interaction of the corresponding layered bulk material. 


3D dielectric function 

ef(q) 


^ D (q) 

q)’ 


(13) 


which is then by construction a diagonal matrix as in¬ 
dicated by the subscript index d labeling the different 
elements on the diagonal. We now aim to link this di¬ 
electric matrix to the interactions taking place between 
electrons within one monolayer and connect it to model 
dielectric functions derived in the context of continuum 
electrostatics. 

To this end, we consider the bare intralayer electronic 
interaction 



z = 0) 


1 




(14) 


and diagonalize it 

^ D (q||) = ■ v?2(q||,z = 0) ■ S(q||), (15) 

which leads to the transformation matrices *5(q||). N qz 
is the number of points used in the (^-summation in Eq. 


m■ The bare intralayer interaction is the same in the 
bulk layered material and in the monolayer, bilayer, etc. 
We assume that the transformation S'(qn) also diagonal¬ 
izes the screened interaction within the 2D system of in¬ 
terest: 


^(q||) = S(q||) • t/J D (q,|) • S f (q||), (16) 


which is indeed the quanity we are searching for. What 
remains is to calculate 


un qn) 


^ D (qn) 

e d D (qn)’ 


(17) 


To obtain we neglect all non-localities of the dielec¬ 
tric screening (c.f. Eq. fl3l) in the vertical direction, i.e. 
we consider e^ D (q||,z = 0) = e^ D (q||, q z )e iq *'° 

and replace the ’’head element” by the model dielectric 
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function from Eq. ED): 


^ D (qii)= (is) 

/ £ eff [Mi D (qib* = °)] \ 

£ i D (qih- = o) 


For simiplicity we also assume 
«d D (qii) = ^ D (qii>* = °) = ^X^ D ( q ii’^) el9 *' 0 - 

qz q z 

(i9) 

Since the transformation matrices <S(q||) and S{ q) are not 
neccessarily the same, Eqs. El and El imply an ad¬ 
ditional approximation, since we neglect this difference 
here. In future works this approximation might be re¬ 
moved by introducing additional transformation matri¬ 
ces. For the graphene systems under consideration, this 
approximation appears to be unproblematic as the com¬ 
parison of cRPA and WFCE data in section HV Dl shows. 

From UV~( qj| ) obtained according to Eq. (fT71) a 
Fourier transformation with respect to qy finally leads 
to screened Coulomb matrices element in real space 



N n 




*qn r ii 


( 20 ) 


which is illustrated in Fig. [3] for the q z = 0 direc¬ 
tion. Here, e is the elementary charge and V is the vol¬ 
ume per atom. The fact that the leading eigenvalue of 
the bare Coulomb interaction in terms of Wannier or¬ 
bitals matches the long wavelength continuum descrip¬ 
tion within large parts of the Brillouin zone very closely 
motivates us to consider exactly this part of the Coulomb 
interaction in the WFCE approach. 

The leading bare interaction eigenvalue is basically in¬ 
dependent of any microscopic properties. This is different 
in the case of the screened interaction. Here, microscopic 
and macroscopic properties are involved through the di¬ 
electric screening of the real material background, as can 
be seen from the right panel of Fig. [31 In the limit of 
small q = |q| —> 0 the tensorial character of the dielec¬ 
tric function becomes obvious. Here, we find en « 3.2 
for the in-plane fields and e±_ ss 2.2 for the out-of-plane 
direction. At larger momentum transfer q , the direction 
dependence of the dielectric function is less pronounced. 
Besides the leading eigenvalues of the Coulomb matrices 
the energetic interval of the other eigenvalues are marked 
by the (light) gray shaded areas in the left panel of Fig. 
[3] for the (bare) screened interaction. These matrix ele¬ 
ments correspond to electronic density variations within 
the unit cell and correspondingly short wavelengths. 


B. Freestanding mono- and bilayer graphene 


which can be used in extended Hubbard models like Eq. 
ED- Here, N Qn is the number of points used in the q^- 
summation. 


IV. FROM GRAPHITE TO GRAPHENE 
HETEROSTRUCTURES 

We derive effective Coulomb interactions of monolayer 
graphene (MLG), bilayer graphene (BLG) and Ir inter¬ 
calated graphite (Gr/Ir) from the interaction calculated 
for bulk graphite in the WFCE approach and benchmark 
our results against direct ab-initio calculations for these 
systems as well as analytical expressions which are valid 
in the long wavelength limit. In all cases, we consider 
Coulomb matrix elements in terms of Wannier functions 
for the carbon p z orbitals. 


We derive the screened Coulomb interactions in free¬ 
standing mono- and bilayer graphene using the bulk 
graphite data in the WFCE approach. The leading eigen¬ 
values of the Coulomb interaction and the correspond¬ 
ing effective dielectric functions are shown as function 
of momentum transfer in Fig. [4j While the compari¬ 
son of results from direct ab-initio calculations and from 
the WFCE approach reveals generally very good agree¬ 
ment, there are some systematic deviations between both 
approaches in the limit of q —> 0. To understand the 
origin of these deviations it is instructive to compare the 
bare Coulomb interaction matrix elements obtained from 
WFCE and ab-intio results to the analytical expression 
for the bare Coulomb interaction between electrons con¬ 
fined to a two-dimensional film in the long-wavelength 
limit: 


A. Bulk Graphite 

The left panel of Fig. [3]shows the leading eigenvalue of 
the bare and screened Coulomb interaction in AB stacked 
graphite, which plays the role of the initial bulk material, 
here. The former is perfectly interpolated by the analytic 
expression 


th(q) 


47re 2 1 

~lT~q 2 


( 21 ) 


V 


2D 

1 


(qy) 


h [+* /h w i , 
2nJ_^ /h V q 2 Qz 

4e 2 arctan (^r) 

A q\\ 


( 22 ) 


where A is the unit cell area per atom and the effec¬ 
tive height h can be chosen to be the interlayer distance 
h = d « 3.35 A for the monolayer and h = 2d for 
the bilayer. At small momentum transfer (quh -C 1), 
this bare interaction approaches the well known limit 
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Figure 3. (Color online) Left: Leading eigenvalues of the bare and screened Coulomb matrix elements of graphite for q z = 0 
obtained from ab-initio calculations together with the analytical description of the unscreened interaction (dashed blue line). 
The (light) gray area indicate the interval of all other eigenvalues of the (bare) screened Coulomb matrix. Right: Momentum 
dependent leading eigenvalue of dielectric function of graphite obtained from cRPA calculations for different values of q z . The 
red markers for q = 0 indicate the parallel (circle) and perpendicular (square) limits of the screening. 
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Figure 4. (Color online) Leading eigenvalue of the dielectric function (outer frame), bare and screened Coulomb matrices (inner 
frame) of monolayer (left) and bilayer graphene (right). Red markers indicate ab-initio calculations and gray markers the 
WFCE values. 


Ui D (q||) —> The term arctan^^-^ in Eq. (l22l) 

plays the role of a “form factor” which accounts for the 
effective height h of the two dimensional layer. For, both, 
the monolayer and the bilayer the WFCE results match 
the analytic expression, which becomes exact for q —>• 0, 
almost perfectly, in contrast to the ab-initio data. The 
ab-initio calculations performed here are in fact supercell 
calculations with periodic boundary conditions. This cre¬ 
ates periodic images of the mono- and bilayer graphene 
in z direction, which contribute to the effective screening. 
We employ an extrapolation to infinite supercell height 
(see Appendix) in order to obtain the monolayer limit, 


which becomes somewhat inaccurate for small q. Thus, 
the deviation of ab-initio and WFCE Coulomb matrix 
elements as well as dielectric functions at small q is likely 
due to this extrapolation problem in the ab-initio data. 

At intermediate q the WFCE and the ab-intio dielectric 
function are in very good agreement. For both, mono- 
and bilayer graphene, the screening rises from 1 to a 
maximum at intermediate q and slightly decreases af¬ 
terwards towards the edges of the first Brillouin zoneii. 
Here, the non-locality (g-dependence) of the screening 
becomes clearly visible. In the long wavelength limit 
the screening vanishes, since we are dealing with a free 
standing two dimensional layer, which is embedded in 










an infinite three dimensional vacuum. By decreasing the 
wavelength, or increasing qy, the Coulomb interaction 
starts to be screened like in a three dimensional bulk 
system, which manifests as an increased value of the di¬ 
electric function. The main differences between the effec¬ 
tive dielectric functions in mono- and bilayer graphene is 
the gradient towards the intermediate maximum and the 
absolute value of the maximum, which are steeper and 
higher, respectively, in the bilayer. I. e. the long range 
Coulomb interaction is less screened in the monolayer 
than in the bilayer, while the short range screening is 
more or less the same. The screened Coulomb interac¬ 
tion obtained from WFCE interpolates the correspond¬ 
ing cRPA data very well, as can be seen from Fig. |4] 
Thus, we have proven, that the WFCE approach to cal¬ 
culate the two dimensional Coulomb repulsion directly 
from the three dimensional bulk data without introduc¬ 
ing additional parameters works very well. 


C. Graphene in a metallic surrounding 
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This system can also be interpreted as Ir intercalated 
graphite. To model this system with the WFCE approach 
we assume perfect metallic screening by Ir, £ 2 , £3 —>• 00 , 
and use the effective height h = 3.35 A of graphene, as 
before. The resulting leading eigenvalues of the Coulomb 
interaction within the carbon p z Wannier orbitals as 
well as the corresponding effective dielectric function are 
shown in Figure [5j The metallic surrounding leads to di¬ 
verging ef D (q) at long wavelengths q —>• 0 in the cRPA 
and in the WFCE approach, as it must be. In contrast to 
free standing mono- and bilayer graphene the screened in¬ 
teractions in Gr/Ir do not diverge at small q, where the 
Coulomb interaction is now efficiently screened by the 
metallic environment. 

The overall characteristics of interactions and screen¬ 
ing as obtained from WFCE agree with the cRPA calcu¬ 
lations. Nevertheless there is a systematic underestima¬ 
tion of the screening ef' D (q) on the order of ~ 17% by 
the WFCE approach as compared to the cRPA. Hence, 
the screened Coulomb interactions are correspondingly 
overestimated by WFCE, here. On physical grounds it 
is clear that the WFCE approach can become inaccurate 
when there is hybridization between e.g. a monolayer of 
graphene and some metallic surrounding. In this case 
the assignment of an effective height h to the graphene 
layer and a separation into a subsystem of graphene and 
’’the environment” is ambiguous. The underestimation 
of the screening in the WFCE approach can indeed be 
cured by decreasing the effective height to h « 2.8 A of 
the modeled monolayer. Treating h as an adjustable pa¬ 
rameter, that is derived from e.g. cRPA calculations is 
one possibility if for instance very complex heterostruc¬ 
tures shall be considered and the intercalated system is 
used as the bulk starting point in WFCE. Here, we are 
taking graphite as the bulk starting point to treat Ir in¬ 
tercalated graphite and keep h = 3.35A to stay with a 
parameter free model. 


Figure 5. (Color online) Leading eigenvalue of the dielectric 
function (outer frame), bare and screened Coulomb matrices 
(inner frame) of Gr/Ir. Red markers indicate cRPA calcula¬ 
tions and gray markers the values from the WFCE approach. 


Regarding the change in electronic interactions, the 
opposite extreme case to going from bulk graphite to free 
standing monolayer is the case of graphene embedded in 
some metallic environment. Perfect metallic screening by 
the environment corresponds to £ 2 , £3 —> 00 , in contrast 
to the case of £2 = £3 = 1 for monolayers surrounded by 
vacuum. In experiment, graphene is frequently grown on 
metals like Ir— or Cu^ or can be surrounded by metals 
e.g. in graphite intercalation compounds^. We consider 
Coulomb interactions in graphene surrounded by Ir in 
the following. 

To this end, we calculated Coulomb interactions for a 
periodically repeated slab composed of graphene mono- 
layer and one “monolayer” of iridium by means of cRPA. 


D. Coulomb Interactions in Real Space 

In order to use the Coulomb terms obtained within 
the WFCE approach in a generalized Hubbard model, 
in which interactions matrix elements enter in real space 
representation, we perform a Fourier transformation: 

< 23 > 

911 uii 

In the case of graphite an additional sum over the q z com¬ 
ponent is performed. The resulting values for density- 
density like £A^(r||) as obtained from cRPA and WFCE 
are given in Table |T| and depicted in Fig. [ 6 ] for mono- 
and bilayer graphene as well as graphite and the Gr/Ir 
system. 

The screened Coulomb interaction in monolayer 
graphene is over the whole ry range bigger than the cor¬ 
responding values of bilayer graphene, graphite, and Ir 
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Figure 6. (Color online) Density-density matrix elements of 
the screened Coulomb interactions for graphite, monolayer 
and bilayer graphene as well as Gr/Ir in real space. Col¬ 
ored markers show ab-initio values, gray markers connected 
by dashed lines show the WFCE values. 


Table I. cRPA and WFCE screened Coulomb interactions for 
graphite, monolayer and bilayer graphene and Ir intercalated 
graphite in eV. Since the AB stacking breaks the sublattice 
symmetry in bilayer graphene for every second neighbor, some 
interactions are given separately for the A and B sublattice 
(in one of the two layers). 


System 

U 0 

Ui 

u 2 

u 3 

Ui 

u 5 

Graphite 

8.1/8.2 

3.6 

2.2/2.2 

1.9 

1.5/1.5 

1.3 

MLG 

cRPA 

9.7 

5.3 

3.8 

3.5 

3.0 

2.8 

WFCE 

9.9 

5.3 

3.8 

3.4 

2.9 

2.6 

BLG 

cRPA 

9.1/9.2 

4.7 

3.2/3.2 

2.9 

2.5/2.5 

2.3 

WCFE 

9.2/9.3 

4.7 

3.3/3.2 

3.0 

2.5/2.5 

2.3 

Gr/Ir 

cRPA 

5.1 

1.5 

0.5 

0.4 

0.2 

0.1 

WCFE 

6.2 

1.8 

0.6 

0.4 

0.2 

0.2 


intercalated graphene. Since the bare Coulomb interac¬ 
tions (not shown here) are nearly the same in all cases, 
variations of the background screened interactions are 
almost entirely due to the successively stronger screen¬ 
ing when going from monolayer graphene via bilayer 
graphene and graphite to graphene encapsulated in a 
metal. 

In agreement with Ref. we find sizeable nonlo¬ 
cal effective Coulomb interactions for graphene, bilayer 
graphene and graphite, which can be however strongly 
reduced due to screening by the environment. This can 
be seen from comparison to the Gr/Ir case. Here, the 
Coulomb interaction is strongly reduced at all r || under 
consideration, i.e. by about a factor of 2 for the local 
terms and more than a factor of 10 for interaction terms 
beyond fourth nearest neighbours. 

The comparison of effective Coulomb interaction ob¬ 


tained from direct cRPA and WFCE calculations shows 
generally very good qualitative agreement with devia¬ 
tions of less than 10%. The only exception are in the 
Gr/Ir data. Here, the local Coulomb interactions are 
overestimated by WFCE by about 1 eV, which is likely 
a result of the approximated effective monolayer height, 
as dicussed in section II V Cl Nevertheless, even in this 
” worst-case” the WFCE approach accounts for ~ 80% of 
the increased screening provided by the metallic environ¬ 
ment. 


V. ELECTRONIC GROUND STATE OF 
BILAYER GRAPHENE HETEROSTRUCTURES 

Bilayer graphene is known to host competing symme¬ 
try broken electronic ground states. Theoretical stud¬ 
ies have predicted that charge- and spin-density waves 
(CDW or SDW), quantum spin hall states (QSH), ne¬ 
matic, superconducting and excitonic insulator states 
could emerge in the bilayer— — . Different experiments 
have addressed the issue of symmetry broken ground 
states in bilayer graphene^ - — but the issue remains con¬ 
troversial also from the experimental point of view and 
it is e.g. unclear whether or not the ground state ex¬ 
hibits a finite electronic excitation gap. In the end, it 
appears very likely that microscopic material specific de¬ 
tails of the effective interactions determine which elec¬ 
tronic phases are realized. 

As a dielectric substrate or also some metallic en¬ 
vironment can provide additional screening of the ef¬ 
fective Coulomb interactions in bilayer, we investigate 
how different types of environments affect the electronic 
ground state of bilayer graphene. To this end we use 
the WFCE approach to study the influence of a dielec¬ 
tric substrate (£3 = 5; £2 = 1) and a metallic substrate 
(£3 —>- 00 ; £2 = 1 ) as well as a metallic encapsulation 
(£2 = £3 —> 00 ). The results can be seen in Fig. [7] as well 
as in Tab. [H] 

The effective dielectric function diverges for q —t 0 in 
the case of the metallic substrate/environment, whereas 
a finite ef D (q = 0) > 1 can be found for the dielectric 
substrate, as it must be. The resulting effective Coulomb 
interactions can be clearly reduced due to environmental 
screening as comparison with the free standing bilayer 
data demonstrates. To understand the resulting effects 
on the electronic ground states, we make use of an elec¬ 
tronic phase diagram that is based on a recent functional 
renormalization group study^ and that gives the ground 
states in dependence of the local, nearest-, and next- 
nearest neighbor screened Coulomb interaction. The in¬ 
teraction strengths obtained here put the ground state 
of bilayer graphene in all environments considered on the 
Antiferromagnetic (AF)-SDW side of a phase transition 
line between the QSH and the AF-SDW phases found 
in Ref. [35. This result of the AF-SDW being stable in 
all situations is quite surprising given the fact, that the 
Coulomb interactions are changed quite drastically. 
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Figure 7. (Color online) Left: WFCE density-density matrix elements of the background screened Coulomb interactions for 
bilayer graphene in real space for different dielectric surroundings (e 2 / £ 3 ). Right: Corresponding dielectric functions in 
momentum space. 


Table II. Screened Coulomb interaction for bilayer graphene 
modified through different dielectric environments £ 2 , 3 - AH 
values are given in units of eV. Since the AB stacking breaks 
the sublattice symmetry for every second neighbor, some in¬ 
teractions are given separately for the A and B sublattice (in 
one of the two layers). 


£2 

£3 

Uo 

Hi 

u 2 

u 3 

1 

1 

9.2/9.3 

4.7 

3.3/3.2 

3.0 

1 

5 

8.2/8.3 

3.7 

2.3/2.3 

2.0 

1 

OO 

7.6/7.7 

3.1 

1.7/1.7 

1.4 

00 

00 

7.3/7.4 

2.9 

1.5/1.5 

1.2 


VI. CONCLUSION 

We have established a scheme that combines cRPA cal¬ 
culations of layered materials in the bulk and continuum 
medium electrostatics to derive effective Coulomb inter¬ 
action matrix elements in terms of Wannier functions 
for free standing 2D materials as well as 2D materials 
embedded in complex dielectric environments. We call 
the scheme ’’Wannier function continuum electrostatics” 
(WFCE) approach. It allows us to avoid supercell calcu¬ 
lations involving complex environments or large vacuum 
volumes on the ab-initio side, which are numerically very 
costly in implementations using periodic boundary con¬ 
ditions. While the WFCE approach might be further 
improved and generalized in the future, already the sim¬ 
plest version presented here predicts effective Coulomb 
matrix elements for monolayer and bilayer graphene very 
accurately, i.e. for instance the local Hubbard interaction 
agrees with the full cRPA calculation within 0.2 eV. 

Our modelling shows that Coulomb interactions can be 
strongly manipulated in 2D materials like graphene by 
means of screening provided by the environments which 
can be substrates, adsorbates or other 2D materials. 


Given the numerical simplicity of the WFCE approach, 
we anticipate that it could be very useful in the context of 
materials design, as effects of different kinds of dielectric 
environments on Coulomb interactions in layered mate¬ 
rials can be modelled quickly and quantitatively, now. 
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Appendix A: Computational Details of the Ab-initio 
Simulations 

In all calculations the carbon nearest neighbor distance 
is set to 1.42 A. In graphite and BLG an AB Bernal stack¬ 
ing is chosen with an interlayer distance of 3.35 AM. To 
model Ir intercalated graphite we use the structure de¬ 
picted in Fig. [HI The distance between the graphene 
layers and the Ir layer is set to the experimental value of 
3.41 AM. 



Figure 8. (a) Top and (b) side view of the Gr/Ir system. 


For all DFT calculations the generalized gradient ap- 
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proximation (PBE)^2 is used. In the case of carbon atoms 
we use an angular momentum cut-off of l cut = 6 and 
lent = 8 for iridium. The plane-wave cut-off is set to 
4.5 a^ 1 , where ao is the Bohr radius. The involved k 
meshes and the energy cut-offs for the polarization func¬ 
tion are shown in Tab. un The energy cut-off corre¬ 
sponds to the energy of the highest, unoccupied band 
(and thus to the total number of empty bands) involved 
in the calculation of the polarization function. Since in 
the case of MLG and BLG several “vacuum distances” 
(see below) have been used, the number of empty bands 
had to be adjusted for each vacuum height (correspond¬ 
ing to the given energy cut-off). 

A “vacuum distance” h vac (the distance between adja¬ 
cent layers) is introduced, since we embed the mono- or 
bilayer in a three dimensional unit cell. Thereby, we pro¬ 
duce, due to the periodic boundaries, an infinite stack 
of mono- or bilayers separated by the unit cell height. 
The freestanding situation is obtained in the limit of 
h v ac —> oo. To approximate this limit, we do several cal¬ 


culations for different vacuum distances (ranging from 
h va .c ~ 15 A to h v ac ~ 30 A) and extrapolate the free¬ 
standing value U a p(q, oo) by fitting the results to 

U a p( q, h v ac) = U a p (q, oo) + , (Al) 


Table III. Ab-initio details for each system. The polarization 
energy cut-offs are given relative to the graphene Dirac-cone 
position. 


system 

DFT k mesh 

cRPA k mesh 

energy cut-off 

AB graphite 

16 x 16 x 5 

25 x 25 x 8 

« 60 eV 

MLG 

14 x 14 x 1 

16 x 16 x 1 

« 120eV 

BLG 

14 x 14 x 1 

28 x 28 x 1 

« 60 eV 

Gr/Ir 

16 x 16 x 5 

28 x 28 x 1 

« 180 eV 
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